Case study

We’ll explore the spatial distribution of ALDI supermarkets in the UK using the sf package. You will need a Mapbox Access Token to use the mapdesk package and generate travel time isochrones.

Load packages

library(sf) # spatial data manipulation 
library(tidyverse) # data manipulation
library(tmaptools) # query OpenStreetMap Nominatim
library(mapboxapi) # access to the Mapbox Isochrone API
library(mapdeck) # interactive maps using mapbox GL and deck.gl

Load data

# Local Authority District boundaries
# Source: ONS Open Geography Portal
# License: Open Government Licence v3.0
# URL: https://geoportal.statistics.gov.uk/datasets/ons::local-authority-districts-december-2021-gb-bgc-1
ltla <- read_sf("data/Local_Authority_Districts_(December_2021)_GB_BGC.geojson")  %>% 
  st_transform(27700) %>% 
  select(AREACD = LAD21CD, AREANM = LAD21NM)

# Lower-layer Super Output Area boundaries
# Source: ONS Open Geography Portal
# License: Open Government Licence v3.0
# URL: https://geoportal.statistics.gov.uk/datasets/ons::lsoa-dec-2021-boundaries-generalised-clipped-ew-bgc-v2
lsoa <- read_sf("data/LSOA_Dec_2021_Boundaries_Generalised_Clipped_EW_BGC_V2.geojson") %>% 
  st_transform(27700) %>% 
  select(AREACD = LSOA21CD, AREANM = LSOA21NM)

# ALDI supermarkets
# Source: Geolytix
# License: CC0 https://creativecommons.org/publicdomain/zero/1.0/
# URL: https://geolytix.com/blog/supermarket-retail-points/
aldi <- read_csv("data/geolytix_retailpoints_v26_202211.csv") %>% 
  st_as_sf(coords = c("long_wgs", "lat_wgs"), crs = 4326) %>%
  st_transform(27700) %>% 
  filter(str_detect(retailer, "Aldi")) %>% # or choose a different retailer
  select(retailer, store_name, town, postcode)

Check data

Get a glimpse of the data

glimpse(ltla)

Rows: 363
Columns: 3
$ AREACD   <chr> "E06000001", "E06000002", "E06000003", "E06000004", "E0600000…
$ AREANM   <chr> "Hartlepool", "Middlesbrough", "Redcar and Cleveland", "Stock…
$ geometry <MULTIPOLYGON [m]> MULTIPOLYGON (((450152.7 52..., MULTIPOLYGON (((…

glimpse(aldi)

Rows: 988
Columns: 5
$ retailer   <chr> "Aldi", "Aldi", "Aldi", "Aldi", "Aldi", "Aldi", "Aldi", "Al…
$ store_name <chr> "Aldi Aberdeen", "Aldi Westhill", "Aldi Ellon", "Aldi Inver…
$ town       <chr> "Aberdeen", "Westhill", "Ellon", "Inverurie", "Hatfield, He…
$ postcode   <chr> "AB11 5EJ", "AB32 6FY", "AB41 9LJ", "AB51 4FY", "AL10 9RQ",…
$ geometry   <POINT [m]> POINT (395153.8 806418.7), POINT (383260.1 807138.5),…

Check that the CRS match

st_crs(ltla) == st_crs(aldi)

[1] TRUE

Exploratory Spatial Data Analysis (ESDA)

Which local authority has the most ALDI supermarkets?

ltla %>% 
  mutate(n = lengths(st_intersects(., aldi)))

How many ALDI supermarkets are there in Greater Manchester?

st_filter(aldi, filter(ltla, AREANM %in% c("Bolton", "Bury", "Manchester", "Oldham", "Rochdale", "Salford", "Stockport", "Tameside", "Trafford", "Wigan")), .predicate = st_within) %>%
  nrow()

[1] 59

Where are the ALDI supermarkets in Cardiff?

cardiff <- filter(ltla, AREANM == "Cardiff")
aldi_cardiff <- st_filter(aldi, cardiff, .predicate = st_within) 

mapdeck(style = mapdeck_style("light")) %>%
  add_polygon(data = st_transform(cardiff, 4326),
              fill_colour = "#00b5db80") %>% 
  add_pointcloud(data = st_transform(aldi_cardiff, 4326),
                 elevation = "z",
                 fill_colour = "#f03d14",
                 tooltip = "store_name") %>% 
  mapdeck_view(location = st_centroid(st_geometry(st_transform(cardiff, 4326)))[[1]], zoom = 11)

How many ALDI supermarkets are within 10 km of the Angel of the North?

point <- geocode_OSM("Angel of the North, England", as.data.frame = TRUE) %>% 
  st_as_sf(coords = c(x = "lon", y = "lat"), crs = 4326) %>% 
  st_transform(27700)

aldi_distance <- aldi %>% 
  st_filter(point, .predicate = st_is_within_distance, dist = 10000) # 10km

mapdeck(style = mapdeck_style("light")) %>%
  add_polygon(data = st_transform(st_buffer(point, dist = 10000), 4326),
              fill_colour = "#00b5db80") %>% 
  add_pointcloud(data = st_transform(aldi_distance, 4326),
                 elevation = "z",
                 fill_colour = "#f03d14",
                 tooltip = "store_name") %>% 
  mapdeck_view(location = st_centroid(st_geometry(st_transform(point, 4326)))[[1]], zoom = 10)

How many ALDI supermarkets are within a 30 minute drive of the Angel of the North?

isochrone <- mb_isochrone("Angel of the North", profile = "driving", time = c(10, 20, 30)) %>%
  st_transform(27700)

aldi_travel_time <- st_filter(aldi, isochrone, .predicate = st_within) 

mapdeck(style = mapdeck_style("light")) %>%
  add_polygon(data = st_transform(isochrone, 4326),
              fill_colour = "time",
              fill_opacity = 0.5,
              legend = TRUE,
              legend_options = list(title = "Time (minutes)")) %>% 
  add_pointcloud(data = st_transform(aldi_travel_time, 4326),
                 elevation = "z",
                 fill_colour = "#f03d14",
                 tooltip = "store_name") %>% 
  mapdeck_view(location = st_centroid(st_geometry(st_transform(isochrone, 4326)))[[1]], zoom = 9)

What is the average drive time to an ALDI supermarket for each of Carlisle’s LSOAs?

carlisle <- filter(ltla, AREANM == "Carlisle")

carlisle_lsoa <- lsoa %>% 
  st_filter(carlisle, .predicate = st_intersects)

aldi_carlisle <- aldi %>%
  st_filter(carlisle, .predicate = st_is_within_distance, dist = 50000)

travel_times <- mb_matrix(carlisle_lsoa, aldi_carlisle)
min_time <- apply(travel_times, 1, min)
carlisle_lsoa$time <- min_time

mapdeck(style = mapdeck_style("light")) %>%
  add_polygon(
    data = st_transform(carlisle_lsoa, 4326), 
    fill_colour = "time",
    fill_opacity = 0.5,
    palette = colourvalues::get_palette("viridis")[256:1,],
    legend = TRUE,
    legend_options = list(title = "Time (minutes)", digits = 0)
  ) %>% 
  add_pointcloud(data = st_transform(aldi_carlisle, 4326),
                 elevation = "z",
                 fill_colour = "#f03d14",
                 tooltip = "store_name") %>% 
  mapdeck_view(location = st_centroid(st_geometry(st_transform(carlisle, 4326)))[[1]], zoom = 8)